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Abstract 

A new scheme for numerical integration of the 1D2V relativistic Vlasov- 
Maxwell system is proposed. Assuming that all particles in a cell of the 
phase space move with the same velocity as that of the particle located 
at the center of the cell at the beginning of each time step, we successfully 
integrate the system with no artificial loss of particles. Furthermore, splitting 
the equations into advection and interaction parts, the method conserves 
the sum of the kinetic energy of particles and the electromagnetic energy. 
Three test problems, the gyration of particles, the Weibel instability, and 
the wakefield acceleration, are solved by using our scheme. We confirm that 
our scheme can reproduce analytical results of the problems. Though we 
deal with the 1D2V relativistic Vlasov-Maxwell system, our method can be 
applied to the 2D3V and 3D3V cases. 

Key words: Relativistic Vlasov-Maxwell system, Plasma instability, 
Laser-plasma interaction 



1. Introduction 

In a wide range of plasma processes operating in laboratories or astrophys- 
ical phenomena, interactions between relativistic particles and electromag- 
netic fields play vital roles. For instance, recent laser experiments revealed 
that a high intensity laser can accelerate particles to ultra-relativistic speed 
(see, e.g., Esarey et al. Non-thermal components found in spectra of 

active astrophysical objects, e.g., supernova remnants, gamma-ray bursts, 
and jets from active galactic nuclei, are interpreted as synchrotron radiation 
emitted by charged, relativistic particles gyrating about magnetic field lines. 
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There are two approaches in order to model such plasmas. One is the fluid 
approach based on the relativistic magnetohydrodynamics(RMHD) and the 
other is the kinetic approach based on the Boltzmann equation coupled with 
the Maxwell equations. Because the fluid approach implicitly assumes that 
the distribution of particles in the momentum space is the Maxwell- Jiittner 
distribution, which is the relativistic extension of the classical Maxwell- 
Boltzmann distribution, the kinetic approach is indispensable for dealing 
with the momentum distribution deviating from the thermal equilibrium. 
Especially, dilute plasmas in which collisions between particles composing 
the plasmas are absent, often called collisionless plasmas, are known to be 
modeled by the so-called Vlasov-Maxwell system [2j. 

At present, the most reliable and reasonable method to simulate dynam- 
ical behaviors of collisionless plasmas is the particle-in-cell (PIC) simulation 
(see, e.g., Birdsall and Langdon which calculates the orbits of charged 
particles by solving the equation of motion and the configuration of elec- 
tromagnetic fields by solving the Maxwell equations. In this method, the 
momentum distribution of plasmas is approximated by an ensemble of the 
momentum of each particle placed in the physical space. Although the num- 
ber of particles in virtual plasmas produced by a PIC simulation is much 
smaller than that in real plasmas, it is known that behaviors of plasmas are 
well reproduced by the method. Nevertheless, we cannot avoid significant 
numerical noises due to the shortage of particles when we focus on the high 
momentum tail of the distribution function of plasma particles. 

On the other hand, the direct numerical integration of the Vlasov-Maxwell 
system (referred to as "Vlasov simulation"), which discretizes the momen- 
tum space as well as the physical space, does not suffer from such noises. 
Therefore, some methods to perform the Vlasov simulation have been de- 
veloped jl, 0, @, 0] . While Vlasov simulations require higher computational 
performance than PIC simulations do, recent developments of computational 
technology allow us to study plasma processes by using Vlasov simulations. 
For example, Mangeney et al. [7] presented a scheme to numerically inte- 
grate the 2D3V Vlasov-Maxwell system (the term "2D3V" means the two 
dimensional space associated with the three dimensional velocity space) and 
demonstrated that the scheme could simulate the Weibel instability with 
high accuracy Valentini et al. [8\ provided a scheme for the integration of the 
electrostatic 1D2V Vlasov-Poisson system in a uniform magnetic field. They 
adopted the polar coordinates in the velocity space, which allows them to 
perform simulations with a good energy conservation. The scheme presented 
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by Valentini et al. 0| integrates the Vlasov- Maxwell system in the hybrid ap- 
proximation, i.e., they solve the 2D3V electromagnetic Vlasov equation for 
ions and fluid equations for electrons, based on the current advance method 



introduced by Matthews [lOfl. Suzuki and Shigeyama [ll| investigated non- 



linear behavior of the Weibel instability in detail by using a scheme similar 



to that of Mangeney et al. [7|. Schmitz and Grauer [12j performed a se- 
ries of simulations for the magnetic reconnection and confirmed that their 
results are consistent with those of some PIC simulations. However, the at- 
tempts stated above treated only non-relativistic plasmas. Investigations into 
the numerical integration of the relativistic Vlasov- Maxwell system are still 
rare. Although Besse et al. [l^] presented a scheme for the 1D2V relativistic 
Vlasov-Maxwell system, they assumed that particles have no dispersion in 
the transverse momentum space. Furthermore, there exists another problem 
that the mass and energy conservations are not always guaranteed unlike 
PIC simulations. 

In this paper, we propose a new conservative scheme for the numerical 
integration of the 1D2V relativistic Vlasov-Maxwell system that allows par- 
ticles to have dispersions in the momentum space. The scheme is based on 
the semi-Lagrangian approach, which is extensively used to solve the Vlasov- 



Maxwell system [13j, [14], uM • m Sec. 121 we introduce the 1D2V relativistic 
Vlasov-Maxwell system and some characteristic scales, and then transform 
the equations for convenience of the subsequent sections. Sec. [3] describes 
the method for the numerical integration of the system. In Sec. IH we calcu- 
late three test problems using the scheme proposed in Sec. [31 the gyration of 
particles, the Weibel instability, and the wakefield acceleration. We conclude 
this paper in Sec. [51 



2. Formulation 

In this section, we present a scheme for the numerical integration of the 
1D2V relativistic Vlasov-Maxwell system. 

2.1. The relativistic Vlasov-Maxwell system 

We consider a plasma whose spatial distribution varies along one direc- 
tion, which implies that only two components of the momenta of particles, 
the longitudinal and the lateral components, need to be calculated. Thus, 
the 1D2V Vlasov equation for species s describes the kinetic evolution of the 
distribution function f s (x,p,q,t) (s = e for electrons and s = i for ions) in 
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the phase space (x,p,q), where x is the coordinate in the physical space, p 
is the corresponding coordinate in the momentum space, and q is the coor- 
dinate labeling the lateral momentum. In this case, the relativistic Vlasov 
equation takes the following form, 



dt m s T s dx \ m s cT s J dp \ m s cT s J dq 

(1) 

where 



2 / \ 2 

V \ , / Q 



1+ — + — (2) 

\m s cj \m s cj 

represents the Lorentz factor. The constants m s and Q s represent the mass 
and the charge of a species s. c is the speed of light. The electric field 
appearing here has two components parallel and normal E 1 - to the 
while the magnetic field has only one component B L normal to the x-axis. 
Here the normal component of the electric field points to the direction of 
the lateral momentum and the electric and magnetic fields are perpendicular 
to each other. Thus, they are expressed as vector forms E = (E^E^^O) 
and B = (0, 0, B- 1 ) when the momentum vector is expressed as p = (p, q, 0). 
Their time evolutions are governed by the Maxwell equations, 

dEW _ _ 47rj|| ldE^ dB^_ _^ jL ldB^ dE^ _ Q 
dt c dt dx c dt dx 

where the electric current densities J" and J 1 - are expressed in terms of 
f s (x,p,q,t) as 

J]l = T,Qs r r -^f'(x,p,q,t)dpdq (4) 

j± = r r - g r-.f a ( x > p > q > t ) dpdq - (s) 

s J -oo J-oo m sCi- 

2.2. Normalization 

For the numerical integration of the equations introduced above, we de- 
fine the characteristic value for each physical quantity: l/w e as the time 
scale, c/uj c as the length scale, m c c as the momentum, cm e u e /e as the elec- 
tromagnetic field, and m e ul/ (4ne) as the electric current density. Here u c is 
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the electron plasma frequency defined by 

47re 2 no 



0J o 



(6) 



where e is the elementary charge and Uq is the number density. Normalizing 
physical variables with these quantities and using the same notations for 
normalized quantities, one can obtain the dimensionless Vlasov equation, 



df s p df s 



EW + ±B- 



dt V s dx 
where the Lorentz factor is modified to 

r 



df s 

+ ^ 



- dq 



Here R s and R s are dimensionless constants defined by 



0, (7) 



(8) 



Rl 



(9) 



respectively. On the other hand, the Maxwell equations lead to the following 
dimensionless form, 



-J 11 



dE 1 - dB J 



dx 



-J J 



dB ± dE 

+ - 



dx 



0. 



(10) 



where the dimensionless electric current densities J" and J 1 - are expressed 
in terms of f s (x,p, q, t) as 



J" 



OO /"OO 



OO J — OO 

oo rca 



^;f s (x,p,q,t)dpdq, 



s 

/oo roo 
/ ^f S (x,p,q,t)dpdq. 
-oo J —oo 



(11) 

(12) 



These two relations close the system. 



2. 3. Transformation of equations 

For convenience of the following sections, we transform equation (I7j) into 
the conservative form and equation ( TTDT) into the advection form. 
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Multiplying equation by T s (p, q) and some algebraic manipulations 
lead to the following equation, 

d(r°n + P d(v s n , Dsf ^ , q dA d(v s n 



+ K ( El1 + fS B ' 



dt T s dx 9 V T s J dp 

P_ B i\ d(T a f a ) _ us pEW+qE^ 
V s ) dq 



+R s E ± _ JL B ±\ -^J—L = Z q f s - (13) 

1 \ VS J Fir, i vs •> v 1 



The L.H.S of this equation represents the advection of the rest and kinetic 
energy of particles along the characteristics of the Vlasov equation (j7]), while 
the R.H.S is interpreted as the exchange of energy between particles and 
electromagnetic fields. 

On the other hand, introducing the following variables, 

G(j) = gM±gjW, g(j) = gCEk-gjM (14) 

one can rewrite the Maxwell equations for the perpendicular components E 1 - 
and B as 

dG dG _J^_ OH__dH_ _J^_ 

~dt + ~dx~ 2' dt dx~ ~2~' ^ ' 

In the following, we integrate the above equations instead of the equations 

for the components E 1 - and B^. 

3. Strategy for numerical integration 

In this section, we describe a method for numerical integration of the 
dimensionless Vlasov-Maxwell system (T7j)- (fT3"l) introduced in the previous 
section. 

3.1. Discretization 

First, we divide the phase space with the range of [x m i n , x max ] x [p mm , p max ] x 
[^min, Q'max] into N x x N p x N q small cells each of which has the volume 
AxApAq. Thus, Ax, Ap, and Aq are 

AX = N x ' AP= iV p ' Ag= iV g • (16) 



6 



The center of a cell labeled by integers (i,j,k) is located at (x,p, q) = 

Pj > 

q k ), where 

a?< = x min + Ax(i - 1/2) for 1 < i < N x , (17) 

Pj Pmin 

+ Ap(j-l/2) for l<j<iV p , (18) 
Qk = <?min + Ag(/c - 1/2) for 1 < k < N q . (19) 

Next, we define the number of particles of a species s in the cell at time t as, 

pXi+Ax/2 ppj+Ap/2 pq k +Aq/2 

N s ijk (t)= dx dp dqf s (x lP ,q,t), (20) 

Jxi-Ax/2 Jpj-Ap/2 Jq k -Aq/2 



and the energy of particles contained in the cell at t; 

PXi+Ax/2 rpj+Ap/2 nq k +Aq/2 

E?jk(t)= dx dp dqrf s (x,p,q,t). (21) 

Jxi-Ax/2 Jpj-Ap/2 Jq k -Aq/2 



On the other hand, we discretize electromagnetic fields by defining them only 
at the positions Xi, 

Ef(t) = E\\( Xi ,t), Ej-(t) = E ± (x u t), B±(t) = B L (x t ,t). (22) 

3.2. Splitting of equations 

Applying the operator splitting method, Equation ffT31) is numerically 
integrated by two steps. One is the step for the advection of particles and 
electromagnetic fields and the other is the step for the exchange of energy 
between particles and electromagnetic fields. 

The Vlasov equation ((7j) is an advection equation with no source term, 
while the energy equation (TT3"|) contains advection terms and a source term. 
We split the energy equation (TT5]) into the two parts as follows, 

d(T°n + p d(T'f) + RS f El{ + q_ B± \ d(T°f 



dt T s dx q \ T s J dp 

+R s ( E ^ _ 1L B ±) ^111 = o, (23) 



dt q V s J K J 

One can see that the advection part of the energy equation fT23|) takes the 
same form as the Vlasov equation (j7|). Therefore, we introduce an operator 
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Ap[E^(t),G(t),H(t),At] that evolves the variables N? jk (t) or E? jk (t) by a 
time interval At according to Equations ([7]) or fl23|) for given E"(t), G(t), and 
(or -E _L (i) and B^it)). For the interaction part, we introduce another 
operator l p [E^(t),G(t), H(t), At] that evolves the variable E^- k (t) by a time 
interval At according to Equation (I21I) with given -E"(t), G(t), and i?(t). 

We present a method to calculate the time evolution of the quantities 
defined by equations (!20|) - (p2|) . As is the case for the energy equation, the 
Maxwell equations contain advection terms and source terms. Thus, we split 
them into the two parts as follows, 

dG dG_ Q ^_^_ (25 ) 

dt dx ' dt dx 

dEW _ || dG dH_ _J^_ 

~dT~ ' dt ~ 2 ' dt ~ 2 ' ( } 

Here we introduce two operators that evolve the variables Gi(t) and Hi(t) 
by a time interval At according to Equations (125]) as ^4 9 [At] and *4.h[At]. In 
addition, for the interaction part, we introduce three operators that evolve 
the variables E-(t), Gi(t), and Hi(t) by a time interval At according to 
Equation (1261) as X e [At], l g [At], and X^[At]. The explicit procedures of the 
thus introduced operators for advection terms are discussed in Sec. 13.31 Sec. 
13.51 discusses those for source terms. Using the operators introduced above, 
we propose a scheme to numerically integrate the relativistic Vlasov-Maxwell 
system according to the following steps, 

stepl : N*„ = A p [EHt),G(t),H(t),At/2]N t s lk (t) 

Et* k = A P lE\\(t),G(t),H(t),At/2}E° 3k (t) 

G* =A g [At]G l (t) 

H* =A h [At]H i (t) (27) 

ste P 2 : E$ = l P [E^t), G*, H*, At]E-* k 

Ef(t + At) =I e [At] 

G? =l g [At]G* 

H** =I h [At]H* (28) 

step3 : N° jk {t + At) = A p [Ef(t + At), G**, H**, Af/2]JV& 
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Ef jk (t + At) = A p [Ej{t + At), GT, H**, M/2]E% 
Gi(t + At) =A g [At)G** 

Hi(t + At) =A h [At]Ht*. (29) 

The electric current densities J 1 - and J", which are necessary for the inte- 
gration of the source terms, are evaluated between stepl and step2. The 
procedure for the evaluation is explained in Sec. 13.41 

3.3. Advection part 

For the integration of the advection part, we make use of the character- 
istics of the Vlasov equation (j7j), 

which are equivalent to the equation of motion of a relativistic charged par- 
ticle, because there exists a reliable scheme for the integration of these equa- 
tions widely used in PIC simulations (if, the Buneman-Boris method. 

At first, using the Buneman-Boris method, we obtain the orbit of a par- 
ticle located at the center of each cell (xi,pj, g&) at time t. We thus calculate 
the coordinates (x'^p'j, q' k ) of the particle at t + At as 

Pj = Pj + 
q'k = Qk + 

We then assume that the other particles in the same cell move with the same 
velocity as the particle having been located at the center, which should be a 
good approximation for a sufficiently small cell. The relation between the size 
of the cell and the accuracy of the above treatment is discussed in §3.7 and 
examined in §4.1. The intuitive explanation for the scheme is shown in Figure 
[TJ In each panel, the horizontal axis represents the x-axis and the vertical 
axis represents the p- and g-axes. Although we draw the phase space as 
two dimensional, actual calculations are performed in the three dimensional 
phase space (x, p, q) . The procedure to calculate the number of particles in 

9 



t+At 



■E-dt, 



t+At 



t+At 



R^[E x -i-B L )dt. 



(31) 



the cell (cell 5) located at the center of the surrounding nine cells at t + At 
is as follows; (1) calculate the orbit of a particle located at the center of each 
cell (the left panel) using the Buneman-Boris method. (2) count the number 
of particles entering the original position of cell 5 under the assumption for a 
uniform distribution of particles inside each cell. In other words, the number 
of particles in cell 5 at t + At is defined as that of particles located in the 
gray zones in the right panel of Figure [TJ Therefore, an explicit expression 
of the operator A p becomes 

i+1 j+1 k+1 

^p[Ef,Gi, Hi, At)N- jk (t) = J2 E E N *M*) 

i>=i-lj'=j-l k'=k-l 

x K ~ x i' \ \Pj'-Pi'\ Wk'-Qk'\ ^ 
Ax Ap Aq 

Here the summations in this expression run over only the cells overlapping 
the original position of the cell at (xi,pj, i.e., cell 5, cell 6, and cell 7 for 
the case of Figured! We evolve the energy contained in a cell E^ k (t) in the 
same way. In this method, the number (or the mass) and the kinetic energy 
of particles are conserved for each step. 

The advection part of the Maxwell equations (1251) consists of two linear 
advection equations with a constant velocity that have exact solutions in the 
form of 

G(x,t)=G(x-t,0), H(x,t) = H(x + t,0). (33) 
Therefore, assuming At = Ax, one finds that the relations 

Gi(t + At) = G f _!(t), Hi(t + At) = H i+1 (t), (34) 

hold. We use these relations for the integration of the advection part of the 
Maxwell equation. Because this method is based on the exact solution of a 
linear advection equation, no numerical diffusion occurs. 

3.4- Interpolation 

As we noted in Section 13.21 the electric current density needs to be eval- 
uated for integration of the interaction part. In the following, we discuss a 
method to evaluate the electric current density. The key ingredient for the 
method is interpolation of the distribution function f s (x,p,q,t). We know 
the number N^ k (t) of particles and the energy Ef- k {t) contained in each cell. 
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From the definition of the two variables, the distribution function f s (x, p, q, t) 
must satisfy Equations ([20]) and (1211) for given Nh k (t) and E-j k (t). In other 
words, we have two constraints. So the interpolation function, which is de- 
fined as fijk{t), generally have the form with two unknown coefficients, 

fijk(t) = a ijk + hjk9(x,p, q), (35) 

where g(x,p, q) is a function and the coefficients a^ k and bijk are determined 
from the constraints (TJOj and (J2TJ) . Here, to determine the form of the func- 
tion g(x,p,q), we consider the meaning of the constraints. The constraints, 
(|20p and f[2Tj) . are the zero-order and the first-order moment of the Lorentz 
factor. Then, we assume the interpolation function flj k (t) to take the form 
of 

f° jk (t)=a ijk + b ijk r s , (36) 

We should note that there are many other candidates for the form of the 
interpolation function. If we calculate time evolutions of other macroscopic 
variable for each cell, e.g., momenta of particles, second-order moment of 
the Lorentz factor, and so on, or use the number N^ ±1 j ±1 ^ k±1 (t) and the 
energy E? ±1 j ±1 k±1 (t) of particles in neighboring cells, we can construct an 
interpolation function including more correction terms, 

fijk{t) = a ijk + kjk? s + c ijk h(x, p, q) H , (37) 

where h(x,p,q) is a function corresponding to the additional macroscopic 
variable. In this study, we use the interpolation function (13"B"1) . which is a 
linear function of the energy of particles, to reduce the computational cost. 

Substitution of the interpolation function (1361) into the constraints and 
some algebraic manipulations lead to 

((n 2 ) Jfc i%(t) - (T s ) Jk Et Jk (t) 



EUt) - (r-) 7 ,.\:>i7) 

AxApAq[((T^} jk -(T%Y 
where the bracket represents the following integral, 



(39) 



i l-Pj+Ap/2 rq k +Aq/2 

( A )i* = "X — ~ / / Mpdq. (40) 



ApAq 
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Appendix IA1 gives the expressions of the variables (T s )j k and ((T s ) 2 )j k . Thus 
the distribution function f?- k takes a uniform value in each spatial cell %. Us- 
ing this interpolation function, the electric current densities due to a particle 
species s are evaluated as 

Pj +Ap/2 r q k +Aq/2 



Jpj-Ap/2 Jq k -Aq/2 L 



RqApAq (a ijk {^) jk + bijktphk) (41) 



Pj+Ap/2 r q k +Aq/2 



fijk = RU / Tr s ft 3 k{t)dpdq 

Jpj-Ap/2 Jq k -Aq/2 L 

= RlApAqiaijki^jk + bijkiq)^) (42) 

in non-dimensional forms. 

However, the thus constructed function f?j k is not guaranteed to take 
positive values at all points in the region \pj — Ap/2,pj + Ap/2] x [q k — 
Aq/2, q k + Aq/2] for each i. Because the distribution function of real plasmas 
must be positive at any point in the phase space, the interpolation is modified 
if f*j k takes a negative value. We use the following simple expressions for jf- k 
and 

- R ^Y~ s) ~Air' ^ ' 

■s± _ ps/ Q \ Njjk(t) (AA s 

3nk - R q\f;/—fa-> ( 44 ) 

instead of the expressions ( )42l) in cells with negative fij k (t). 

One can evaluate the electric current density by summing up these vari- 
ables as 

4 = EEE>& ^EEE* («) 



3. 5. Interaction part 

In this subsection, we propose a method to integrate the interaction part 
with respect to time. This method conserves the sum of the kinetic energy 
of particles and the electromagnetic energy. 
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Equations are discretized as 



E}(t + At) 
Gi(t + At) 

Hi(t + At) 



Ef(t) - j\t + At/2) At 



At, 



(46) 



where the electric current densities are evaluated beforehand according to 
the procedure described in the previous subsection. These equations give 
expressions for the operators Z e , X g , and 2^. On the other hand, to obtain 
the energy E?j k (t) of particles in a cell evolved by the interaction part of the 
energy equation (1241) . Equation ( 1241) is integrated with respect to p and q as 



5_ 

dt 



Pj+Ap/2 fq k +Aq/2 



Pl-Ap/2 Jq k -Aq/2 



T s f s dpdq 
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Pj+Ap/2 r q k +Aq/2 



Pj-Ap/2 Jq k -Aq/2 



P 



f s dpdq (47) 



+E 



± 



Pj+Ap/2 r q k +Aq/2 



Pj-Ap/2 Jq k -Aq/2 



— f s dpdq, 



which means that the total energy of particles in a cell is changed by inter- 
actions between particles and electric fields. Discretizing this equation, we 
then propose the following scheme for the integration: 

E s ijk {t + At) - E s ijk {t) E\t + At) + £11 (t) . s „ , E±(t + At) + E^t) 



At 



-i S± 

Jijki 

(48) 

which gives an expression for the operator X p . 

In the following, we will show that this procedure conserves the total 
energy. Summing up the above equation with respect to j, k, and s, and 
then substituting the relations (113]) and (1151) into the result, one obtains 



jks jks 



E }(t)-^P-At 



+ 



Gi^+Hi^-^-AlAt 



4(t)At (49) 
Jt{t)At, 
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which represents the change of kinetic energy of particles after a time step 
in this scheme. The change of the electromagnetic energy is obtained by 
summing the square of each of (T4*6l) . 



2 

. [^!(*)] 2 



[G t (t)} 2 + [HiW 



T±(t) 

Gi(t)+Hi(t)--^At 



J*(t)At 
Jf(t)At (50) 



By summing up both sides of Equations ( 1491) and ( 1501) with respect to i, one 
can easily check that the total energy in a region [xi — Ax/2,Xj + Ax/2] x 
[PmimPmax] x [g min , g max ] is conserved: 

+ + [GM 2 + mt)f = const. (51) 

ijks 

In other words, the procedures T e , I g , X^, and X p expressed by ( )46|) and ( )48|) 
give a conservative scheme for the integration of the interaction part of the 
relativistic Vlasov-Maxwell system. 

3. 6. Conditions for the time interval 

In Sections 13.21 13. 3[ 13.41 and 13. 5[ we present procedures that evolve the 
number and the energy of particles in a cell and electromagnetic fields. In 
order for the procedures to work, the time interval At is required to satisfy 
some conditions. 

As we noted in Section [3TB| the scheme (I3~4"]) requires that the time interval 
At must be equal to Ax. Furthermore, the scheme for integration of the 
advection part of the Vlasov equation (132]) requires that the displacement of 
a particle by integration of Equations (!30|) along the x-, p-, and g-axes must 
not exceed the intervals Ax, Ap, and Ag. In short, particles must not jump 
over a cell. These conditions impose the value of the time interval to satisfy 

At = Ax, At < ^ , At < Aq , — — , (52) 

max(|^J| + |^|)' max(|^| + |^|)' 1 ; 

where max(Aj) represents the maximum of the variable Ai for all i. 
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3. 7. Accuracy of the scheme 

Finally, we mention the accuracy of our scheme proposed in this section. 
As explained above, our scheme is based on various procedures, such as split- 
ting of equations, the advection part, the interaction part, and the evaluation 
of the current density, which makes the mathematical proof of the accuracy 
of our scheme very difficult. Then, we estimate the accuracy of the advection 
of particles, which is likely to be the most inaccurate compared to the other 
procedure. 

In the procedure solving the advection part of the Vlasov equation, all 
particles in a given cell in the phase space are assumed to move with the 
same orbit as that of the particle located at the center of the cell. However, 
this treatment obviously involves errors to a certain extent, because particles 
located at the different position from the center must be integrated under 
different initial conditions. In particular, the difference is most significant for 
particles located at the vertex of the cell. Since the difference of the position 
between particles at the vertex and the center is of the order of Ax, Ap, and 
Aq, the estimated positions of particles at the vertex contain errors of the 
order of Ax, Ap, and Aq, which indicates the number of particles in the cell 
at the next step contains errors of the order of Ax, Ap, and Aq. Therefore, 
the procedure solving the advection part of the Vlasov equation have first 
order accuracy in the physical and the momentum spaces. 

4. Test problems 

In this section, we show results of simulations performed by using the 
scheme proposed in the previous section. For the purpose, we solve three 
test problems, the gyration of particles, the Weibel instability, and the wake- 
field acceleration. The gyration of particles is solved to investigate into the 
accuracy of our scheme. The Weibel instability and the wakefield accelera- 
tion are well-known plasma processes and important in both experimental 
and astrophysical contexts. 

4-1- Gyration of particles 

We assume that electrons are uniformly distributed in the physical space 
with a gaussian distribution in the momentum spaces, 




(53) 
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where a represents the dispersion in the momentum spaces, and an uniform 
magnetic field, 



where B is a constant. One can easily check that the above configuration 
is a stationary solution of the Vlasov-Maxwell system. However, since our 
scheme suffers from a numerical diffusion as expected in the previous section, 
the distribution function f e (x,p,q,t) at t must be different slightly from the 
initial one f e (x,p,q,0). Then, we adopt the following value e as a measure 
of the accuracy of our scheme, 



where t g represent the gyration period given by m e c/(eBo). Figure [2] shows 
the result with a 2 = 2.0 and Bq = 1.0. The ranges of the space coor- 
dinate, the longitudinal momentum, and the lateral momentum are given 
by x G [— V2tt, v27r], P, q € [—10,10], respectively. The periodic bound- 
ary is imposed in the physical space, while, in the momentum space, the 
free boundary condition is imposed. The filled circles represent the values e 
for various N x (= 110, 120, 130, 140, 150, 160, 170, 180, and 190) and fixed N p 
and N q (N p = N q = 200), whereas the filled squares represent those for 
N p = 110, 120, 130, 140, 150, 160, 170, 180, and 190 and N x = N q = 200. 

The solid line shows that the value e seems to scale roughly as (Ax) L8 . 
The value e is expected to depend strongly on Ap and Aq rather than Ax 
in this test problem where particles rotate in the momentum space (p,q). 
In other words, the dependence of e on Ax have uncertainty because of the 
insensitiveness. Therefore, we conclude that the dependence derived above 
must be e oc (Ax) 2 essentially. However, this does not mean second order 
accuracy in the physical space. Because of the condition At = Ax mentioned 
in §3.6, when we double the number of zones N x in the physical space, the 
time interval At must be half of the previous value. Therefore, the value e 
scales as AtAx, which indicates that our scheme has first order accuracy in 
time and the physical space. On the other hand, the dashed line shows that 
the value e scales as Ap, which confirms the estimation in §3.7. 

4-2. Weibel instability 

The Weibel instability is a kind of plasma instabilities caused by anisotropic 
momentum distributions of collisionless plasma. The formulation and dis- 



£%,0) = E ± (x,0) = 0, B ± (x,0) = B t 



(54) 




(55) 
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persion relation of the Weibel instability operating in a relativistic one- 
dimensional plasma are shown in Appendix [Bj 

For a simulation of the Weibel instability, we treat ions as a uniform 
background and assume that electrons have the following initial distribution, 

,%-9b) +5(q + q h ) 



f e (x,p,q,0) = 5(p)- 

which is approximated in the discretized form by 
1/2 for 



(56) 



1/2 



for 



otherwise 



for 



v/TT^/2 for 







otherwise 



Ap/2 < pj < Ap/2 

and q h - Aq/2 < q k < q h + Aq/2, 
Ap/2 < pj < Ap/2 
and - q h - Aq/2 < q k < -q h + Aq/2, 



(57) 

-Ap/2 < pj < Ap/2 

and q h - Aq/2 < q k < q h + Aq/2, 
-Ap/2 < pj < Ap/2 

and - q h - Aq/2 < q k < -q h + Aq/2, 

(58) 



Here we have introduced a constant % that represents the bulk momentum 
of the counter-stream of the plasma. The initial configuration of the electro- 
magnetic field is 

E^=E ± = 0, B L = ecos(kx), (59) 

where e is a small parameter (= 10~ 5 ) and k is the wave number of the 
perturbation. 

We calculate the evolution of a plasma with the above initial condition 
in the simulation domain whose spatial interval is given by x G [—7i/k,7i/k]. 
The longitudinal momentum ranges are given by p G [—5, 5] for % — 2.065 
(the corresponding bulk velocity is 0.9c), p G [10, —10] for q h = 7.018 (0.99c), 
and p G [30, —30] for % — 22.344 (0.999c). The lateral momentum range is 
q G [—5,5]. The periodic boundary is imposed in the x direction, while, in 
the momentum space, the free boundary condition is imposed. 

Figure [3] shows the time evolutions of the kinetic energy K e of electrons, 
the electric energy E, and the magnetic energy B defined by 



K P 



ijk 



(*)> 



(60) 
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i 

= ^EN') 2 +(^+^) 2 ]> ( 61 ) 

i 

B = ^E(^) 2 = ^E(^-^) 2 ' ( 62 ) 

i i 

for the case of q^ = 2.065 (the corresponding bulk velocity is 0.9c) and 
k — 1. The numbers of zones for the three coordinates are N x = 100, 
N p = N q = 50 for q h = 2.065, JV a = 100, N p = 100, N q = 50 for g b = 7.018, 
and N x = 100, N p = 300, JV 5 = 50 for q h = 22.344. The dashed line in 
Figure [3] reproduces the growth rate calculated from the linearized analysis 
with P t h = 0.1 described in Appendix [B] Although we treat a cold plasma 
whose initial momentum distribution is given by (|56p . the initial setup (|57|) 
has particles with some dispersions in the momenta of the order of the width 
of the momentum bins Ap = Aq = 0.2. Therefore we compare the dispersion 
relation from the numerical simulations with that derived from linearized 
analyses with a finite temperature corresponding to the size of the momentum 
bin. The numerical simulation seems to reproduce the theoretical growth rate 
for a given wave number k — 1. 

Figure H] summarizes the growth rates for other cases as a function of 
the wave number of the perturbation. The lines in this figure represent the 
growth rates for the bulk velocities 0.9c, 0.99c, and 0.999c calculated from 
the dispersion relation (19~4"]) with P th = 0.1. The plotted points show values 
measured from results of the simulation. 



4-3. Wakefield acceleration 

The wakefield acceleration is a promising mechanism for the acceleration 
of particles to highly relativistic speeds (see, Esarey et al. [lj, for review). 
The ponderomotive force of a coherent electromagnetic wave propagating 
in a stationary plasma, such as an intense laser in laboratory or a light 
pulse emitted by a certain active phenomenon in astrophysical environment 



[161 . 1171 . |18| , excites a longitudinal electric field and efficiently generates high- 
energy particles. 

To simulate such circumstances, we impose the following boundary con- 
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dition for the electromagnetic field, 



G(0, t) = A u>L exp 



(t - 2i 



T 2 



sin(o; L t), H(0,t)=0, (63) 



which produces a linearly polarized electromagnetic wave (light pulse) prop- 
agating in the +x-direction. Here we have introduced some parameters char- 
acterizing simulations; A the scale of the vector potential, the frequency 
of the light pulse, r the duration of the light pulse. For the initial con- 
figuration of particles, we consider a cold, homogeneous, stationary plasma 
composed of electrons. The momentum distribution is expressed as 



f e (x,p,q,0)=6(p)6(q), (64) 



which leads to 



1 for -Ap/2 < pj < Ap/2 
N° jk (0) = Et jk (0) ={ and - Aq/2 < q k < Aq/2, (65) 

otherwise 

We treat ions as a neutralizing background, choosing the value of the scale of 
the light pulse to be longer than the electron inertial length c/ u e but shorter 
than the ion inertial length m{/ m & c/ uj e . 

Some results for the case of A = 2.0, = 2.0, and t = n/2 are shown 
in Figures [5], EJ and [3 In each figure, the panels represent the color-coded q- 
integrated distribution function, the longitudinal electric field, the transverse 
electric field, and the transverse magnetic field from top to bottom. It is 
clearly seen that a sinusoidal electrostatic field is excited immediately after 
the passage of the light pulse and accelerate electrons, resulting in some 
bunches of electrons in the phase space. The responses of the plasma and 
the electric field to the light pulse is consistent with the previous studies. 
Sprangle et al. jio| studied this process by numerically solving equations 
which treat non-linear interactions of particles and waves (see also, Ting et 
al. [21f| ) . They showed that sawtooth-like longitudinal waves associated with 
some bunches of particles in the phase space form after the passage of a 
light pulse. Recent two-dimensional PIC simulations (see, e.g., Kuramitsu 



et al. [19|) show a similar behavior. The behavior is also reproduced in our 
results, which indicates that the method presented here can treat the correct 
behavior of the distribution function of relativistic, collisionless plasmas. 



19 



The electron distributions in the momentum space at t = 200 and x = 
178, 184 are plotted in Figure [HI It is clearly seen that the existence of high- 
energy electrons up to p = Ylm c c (the corresponding velocity is equal to 
0.998c) at x = 178, where the strong electrostatic field, i.e., the wakefield, is 
excited due to the ponderomotive force. At x — 184, on the contrary, there 
exists no accelerated electron, since the node of the wakefield is located at 
this point. 



5. Discussion and Conclusions 

In this paper, we have proposed a new conservative scheme for numerical 
integration of the relativistic Vlasov-Maxwell system and performed three 
test problems, the gyration of particles, the Weibel instability, and the wake- 
field acceleration. Adopting semi-Lagrange method, we succeed in developing 
a scheme that conserves the number of particles and the sum of the ener- 
gies of particles and electromagnetic fields. Since the previous scheme [3] 
solving the relativistic Vlasov-Maxwell system do not treat the dispersion of 
the lateral momentum of particles, our scheme is the first one that can treat 
the dispersion correctly. Results of the simulations clearly indicate that our 
method succeeds in reproducing detailed behaviors of the distribution func- 
tions in the phase space. Especially, the tail of the distributions where only a 
tiny fraction of particles reside seems to be solved with considerably high ac- 
curacy, while PIC simulations would suffer from large statistical error there. 

As we noted above, Vlasov simulations generally require more computa- 
tional resources than PIC simulations do. Furthermore, as previous works 
@, IS] have investigated, Vlasov simulations suffer from so-called " filamen- 
tation problem" . Ref. [23J studied wave-particle interactions of a plasma ap- 
proaching to an equilibrium state using PIC simulation and showed that the 
equilibrium is realized through a phase mixing accompanied by formation of 
filamentary structures in the phase space. In Vlasov description, particles 
composing a plasma are treated as a continuous medium, which means that 
Vlasov equation can not take account of essential discreteness of plasma. As 
a result, artificial entropy may arise when a structure with the characteristic 



scale smaller than the mesh size is generated in the phase space. Ref. [23 
argued that this artificial entropy prevents the plasma treated by the Vlasov 
simulation from following the correct path toward the statistical equilibrium. 
Therefore, we need careful studies of long-term evolutions of plasmas by using 
Vlasov simulations. 
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Nevertheless, Vlasov simulations provide us detailed dynamics of plasmas 
in the phase space. Though we deal with the 1D2V relativistic Vlasov- 
Maxwell system, our method can be applied to the 2D3V and 3D3V cases 
Although our scheme proposed in this paper suffer from numerical diffusion, 
there is a plenty room for improvement. For example, in order to integrate 
advection part of the Vlasov equation, we can use the orbit of the particle 
located at each vertex of a cell. In other words, taking account for the 
deformation of the cell at each time step improves the accuracy of the scheme. 
In theoretical investigations into complex behaviors of collisionless plasmas, 
Vlasov simulations must be a attractive tool to compensate defects of PIC 
simulations. 
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A. Evaluation of Some Integrals 



In this section, we evaluate some integrals used to construct the interpo- 
lation function in Sec. 3. 4. The first one is the average of the modified Lorentz 
factor (jSj) in the phase space, (T s )jk, defined by 

i fPj+Ap/2 pq k +Aq/2 

(Hi* = 7 ~ W - / / V(R s J 2 +P 2 + q 2 dpdq. (66) 

LAplAq J pj -& p /2 Jq k -Aq/2 

To perform the integrations, we introduce a function of the variables p and 
q described by 



f V(R s J 2 +p 2 + i 2 



t^Arctan 



pq 



R° m V(R s J 2 + p 2 + q 2 



+|[p 2 + 3(C) 2 ] In [q + VW+pH? 2 



+|[g 2 + 3(0 2 ]ln p+ ^(RsJ2 +p 2 + q 2 
6 L 



(67) 



Since the differentiation with respect to p and the subsequent differentiation 
with respect to q of this function leads to 

d 2 9l 



dpdq 



^(R s m ) 2 +p 2 + q 2 



(68) 
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one can evaluate the integral (r s ) as 

9l ( Pj + Ap/2, q k + Aq/2) 9l ( Pj + Ap/2, q k - Aq/2) 



(n 3k = 



ApAq ApAq 
g 1 (p j - Ap/2, q k + Aq/2) + gi ( Pj - Ap/2, q k - Aq/2) 



ApAq ApAq 

The second integral is the average of the square of the modified Lorentz factor 
in the phase space, ((T s ) 2 ) jk , defined by 

i rpj+Ap/2 rq k +Aq/2 

« P ) 2 ^ = / / [i RS m) 2 +P 2 + 1"] dpdq- (70) 

l±Pl±q J p,- Ap/2 Jq k -Aq/2 

This integration is straightforward and one obtains 



<(n 2 >, fc = (/4) 2 + P 2 + ^ + £ + ( 71 ) 

The third and forth integrals are defined by 

, „ \ 1 /"Pj+Ap/2 pq k +Aq/2 

(JL) = -±— / / P =dpdq, (72) 

and 

g , ]_ r Pj +Ap/2 rqk+Aq/2 



rq k -t£±q/z 

/ rfprfg, (73) 

Jq k -Aq/2 \ (EL) 2 + V 2 + O 2 



X*/ jfe ApAqJ Pj _ Ap/2 J qk - Aq/2 ^J(R s m ) 2 +p 2 + q 2 
respectively. To evaluate the third integral, we define the following function, 

92(p, q) = lV(R s J 2 +P 2 + Q 2 + {RSJ l + ^ In [q + V + P 2 + q 2 ■ 

(74) 

Since the differentiation with respect to q and the subsequent differentiation 
with respect to p leads to 

(dg 2 \ p (?5) 



dp\dq) v /(i^)2 +p 2 + ? 2' 

the third integral is written as 

p\ g 2 (p 3 + Ap/2, q k + Aq/2) g 2 ( Pj + Ap/2, q k - Aq/2) 



V s I jk ApAq ApAq 

g 2 {p 3 - Ap/2, q k + Aq/2) + g 2 ( Pj - Ap/2, q k - Ag/2) 



ApAg ApAg 
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Using the same function, the forth integral is expressed as 

q V g 2 {q k + Aq/2, Pj + Ap/2) g 2 (q k - Aq/2, Pj + Ap/2) 



V s I jk ApAq ApAq 

g 2 (q k + Aq/2, Pj - Ap/2) + g 2 (q k - Aq/2, Pj - Ap/2) 



ApAq ApAq 

The remaining integrals (p)j k and (q)j k can be evaluated by a straightforward 
manner: 

(P)jk = Pj, (q)jk = qk- (78) 

B. The Relativistic Weibel Instability 

The dispersion relation of the Weibel instability in both nonrelativistic 
and relativistic plasmas have already been derived in several investigations 
(see, e.g., Califano et al. |22|). Nevertheless, we review the formulation 
and the dispersion relation of the Weibel instability in a relativistic one- 
dimensional plasma for completeness of this paper. 

B.l. Formulation 

We assume that the initial state characterized by an unperturbed distri- 
bution function /q has no electromagnetic field and that the space (x) and 
time (t) dependences of the perturbations are proportional to exp[i(kx — uit)). 
We then consider how the perturbation 5f s on the distribution function and 
the lateral components of the electromagnetic fields 5E ± and 5B ± evolve ac- 
cording to the relativistic Vlasov-Maxwell system. The linearized relativistic 
Vlasov equation expressed as 

(-fc, + »£) 6f + K^Jl + R- (SE- - f t SB-) M = 0, (79) 

and the linearized Maxwell equations 

-iujSE 1 - + ik5B L = -SJ^, 

-iuj5B L + ik5E ± = (80) 

govern the time evolutions of the perturbed quantities. Elimination of 5B 1 - 
in Equations ( 1791) by using Equation ( IHUl) yields 

V T s / J q T s uj dp q \ V s k J dq v ; 
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and the expression for Sf s , 



iR l ( * dfg , dfS\ SI , ± 



The perturbed electric current density 5 J 1 - is related to the perturbed dis- 
tribution function 5f s as 

5 J 1 - = R S q [°° dp r dq^Sf, (83) 

J— OO J —CO 

and necessary to obtain the dispersion relation. We then assume that ions 
are uniformly distributed in the physical space with no bulk velocity, 

ti = n 6(p)6(q), (84) 

where S(x) represents the delta function, and that electrons have the following 
form of the initial distribution, 

= no e(p + P lh )-e ( p-P tt ) 

Q(q + P + P th ) - 8(g + P - P th ) + Q(q - P + P th ) - 8(g - P - P th ) 

4P th 

where 0(x) represents the Heaviside function. The parameter P t h represents 
the thermal dispersion of the momentum distribution of electrons, and Po 
represents their bulk momentum. This assumption means that only electrons 
contribute to the generation of the electric current density. We define the 
following two integrals, 

/, = r dp r ^, (86) 

which contribute to the electric current density 5J ± , and evaluate them be- 
forehand. 

From the properties of the Heaviside function, the first integral reduces 

to 

n k 2 f Po+p ^ q 2 dq 



h = 
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We then define the function g(q) in the form of 



9(0) 



uj 2 (1 + Pi) - PP t 2 h Arctan 



kR 



th 



-kP th log(q+Jl + Pl + q 2 



Since the derivative of this function is 



u 2 kP 



th 



dg _ 

dq " v /l + ^ h + ? 2 ^( 1 + P l + <? 2 ) - ^th' 



y/u>*(l + Pl)-k*Ply/l + Pl + q* 

(89) 

(90) 



one can express the integral 1\ as 

n k 



h 



2*> 2 



[g(P + P th )-g(P -P th )}. 



(91) 



On the other hand, the second integral becomes 



n 



2P 2 



[P + P th )Arcsinh 



th 



y/1 + (P + ^th) 



-(P - P th )Arcsinh 



P 



th 



1 + (P - P th ) 2 



(92) 



from the straightforward evaluation of the integral. 

Using these integrals, the perturbed electric current density is expressed 

as 



5J J 



10 



ih + I 2 )8E ± . 



(93) 



Substitution of this expression into Equation flHUj) and non-trivial 5B 1 - yields 
the following dispersion relation, 



u 2 -k 2 + ^(I 1 + I 2 ) = 
n 



(94) 



B.2. Dispersion relation 
B.2.1. cold plasmas 

Before we solve the dispersion relation ( |94l) with a fixed wave number k 
and obtain the frequency u, we simplify Equation ( |94|) by taking cold limit 
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(-fth ~~ > 0) to clarify whether any unstable mode exists or not. Under the 
cold limit, the integrals defined in the previous section becomes 



and 



lim Ii 



lim J 2 



n k 2 



p2 
M) 



;i + p 2 )3/2- 



n 



a+p 2 ) 3/r 



respectively. Then, the dispersion relation becomes 

1 



k 2 + 



i + p 2 ) 3 / 2 



u 2 - 



k 2 P 2 



(1 + P 2 )3/2 



0. 



(95) 



(96) 



(97) 



One of the solution of this equation is a pure imaginary number, which means 
that this dispersion relation contains at least one unstable mode. The growth 
rates of the mode 7 defined by 27 = u versus wave number k for initial bulk 
velocities 0.9c, 0.99c, and 0.999c are plotted as solid lines in Figures [H [TOj 
anddU 



B.2.2. warm plasmas 

The analysis performed in the previous subsection indicates that there 
exists an unstable mode. Solving the dispersion relation of warm plasmas 
(I94p . one obtains the growth rate of the Weibel instability for plasmas with 
finite temperature. The results for the case of P t h = 0.1 and the bulk ve- 
locities 0.9c, 0.99c, and 0.999c are also shown in Figures [91 [101 and fTTT One 
can see that the growth of unstable modes is suppressed in the high wave 
number regime. 
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Figure 1: Schematic views of the integration of the advection part. 
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Figure 3: The time evolution of the electron energy, K e (solid), the electric energy E (dash- 
dotted), and the magnetic energy B (dashed). The dashed line represents the theoretical 
growth rate derived by the linearized analysis. 
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Figure 4: The dispersion relation of the Weibel instability. The horizontal axis represents 
the wave number of perturbation and the vertical axis represents the corresponding growth 
rate. The solid, dashed, and dotted lines corresponds to the case that the bulk velocity is 
0.9c, 0.99c, and 0.999c. The points plotted the plane is the value measured from results 
of the simulation. 
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Figure 5: Snapshot of the distribution function and the electromagnetic fields at t = 10. 
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Figure 7: Same as Figure [5j but for t = 200. 
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Figure 8: The g-integrated electron distribution in the longitudinal momentum space at 
t = 200 and x = 178, 184. 
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Figure 9: The growth rate of the relativistic Weibel instability with the bulk velocity 0.9c 
as functions of wave numbers. The solid line represents the growth rate calculated from 
the dispersion relation for cold plasmas ([97]). The dashed line represents the growth rate 
calculated from the dispersion relation for warm plasmas with P t h = 0.1 ([94)) . 
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Figure 10: Same as Figure [9l but for the bulk velocity 0.99c. 
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Figure 11: Same as Figure [9l but for the bulk velocity 0.999c. 
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